!************************************************************************************************************************************
! Much of the code for mIntegration is based on Numerical Recipes routines,
! but it is arranged in a more object-oriented fashion.
! The default is to perform a Romberg integration. There are some functions for quadrature methods, but I haven't used
! them much and they should probably be tested.
!
! The way to use it is to declare an integral object:
!
! type(tIntegral) :: integral
!
! and call the constructor:
!
! new(integral, ll, ul, [eps], [open], [changeVariables], [params], [integerParams],
!               [char], [universe], [lps], [halo], [smith], [galaxyBin], [chi])
!
! You must specify the integration limits: ll, ul.
!
! The other arguments are optional:
!
! [eps] sets the (relative) tolerance to declare convergence.
! [open] is logical and picks open romberg integration (doesn't not evaluate the function at the integration limits)
! [changeVariables]: 'EXP' is most useful. If you would plot the integrand on a log scale, try 'EXP'.
!                    Also, 'INV' is good when one of the limits is infinity.
! [params], [integerParams], [char], [universe], [lps], [halo], [smith], [galaxyBin], [chi]:
!      All these exist to allow values or types to be passed into the integrand function. You can access them with get methods, like
!     lps = getLPS(integral)
!
! Once constructed, perform the integration by calling the function as
!
! call integrate(integral, integrand, answer)
!
! where integrand is the function you want to integrate. This function must have the interface integrand(integral, x),
! where x is a real(D).
!
! You can then check that the integral converged by calling:
!
! call checkConverged(integral, errorMessage)
!
! If it didn't converge to eps tolerance, the program will stop with the message errorMessage (character).
! There are good examples in base/linearPower.F90.
!
!************************************************************************************************************************************
module mIntegration

        use mPrecision
        use mObject

        implicit none

        interface new
                module procedure new_integral
        end interface

        interface integrate
                module procedure integrate_0
                module procedure integrate_2
        end interface

        interface getUniverse
                module procedure getUniverse_integral
        end interface getUniverse

        interface getParam
                module procedure getParam_integral
        end interface getParam

        interface getParams
                module procedure getParams_integral
                module procedure getParams_2
                module procedure getParams_3
                module procedure getParams_5
                module procedure getIntegerParams_2
                module procedure getIntegerParams_3
        end interface getParams

        interface getLPS
                module procedure getLPS_integral
        end interface getLPS

        interface getSmith
                module procedure getSmith_integral
        end interface getSmith

        interface getGalaxyBin
                module procedure getGalaxyBin_integral
        end interface getGalaxyBin

        interface getChi
                module procedure getChi_integral
        end interface getChi

        interface converged
                module procedure converged_integral
        end interface converged

        interface checkConvergence
                module procedure checkConvergence_integral
        end interface checkConvergence

        private
        public :: new
!        public :: converged, checkConvergence, converged_2
        public :: checkConvergence, converged_2
        public :: integrate, quadrature, gaussLaguerre
        public :: getParam, getParams, getChar, getUniverse, getLPS, getHalo, getSmith, &
                  getChi, getGalaxyBin

contains

        subroutine new_integral(integral, ll, ul, eps, open, &
            changeVariables, params, integerParams, char, universe, lps, halo, smith, &
            galaxyBin, chi)

                type(tIntegral), intent(out) :: integral
                real(D), intent(in) :: ll, ul
                real(D), intent(in), optional :: eps
                logical, intent(in), optional :: open
                character(len = 3), intent(in), optional :: changeVariables
                real(D), intent(in), optional :: params(:)
                integer, intent(in), optional :: integerParams(:)

                real(D), parameter :: INTEGRAL_DEFAULT_EPS = epsilon(1.0_R4)

                character(len = INTEGRAL_CHAR_LEN), intent(in), optional :: char
                type(tUniverse), intent(in), optional :: universe
                type(tLPS), intent(in), optional :: lps
                type(tHalo), intent(in), optional :: halo
                type(tSmith), intent(in), optional :: smith
                type(tGalaxyBin), intent(in), optional :: galaxyBin
                type(tChi), intent(in), optional :: chi

                integral%ll = ll
                integral%ul = ul

                if( present(eps) ) then
                        integral%eps = eps
                else
                        integral%eps = INTEGRAL_DEFAULT_EPS
                end if

                integral%dimension = 0

                integral%converged = .false.
                integral%converged_2 = .false.
                integral%error = .false.

                if( present(open) ) then
                        integral%open = open
                else
                        integral%open = .false.
                end if

                if( present(changeVariables) ) then
                        integral%changeVariables = changeVariables
                else
                        integral%changeVariables = ''
                end if

                if( present(params) ) then
                        if( size(params) > INTEGRAL_MAX_PARAMS ) then
                                write(0, *) 'ERROR: Too many parameters given to integral'
                                stop
                        else
                                integral%params( 1:size(params) ) = params
                        end if
                end if

                if( present(integerParams) ) then
                        if( size(integerParams) > INTEGRAL_MAX_INTEGER_PARAMS ) then
                                write(0, *) 'ERROR: Too many integer parameters given to integral'
                                stop
                        else
                                integral%integerParams( 1:size(integerParams) ) = integerParams
                        end if
                end if

                if( present(char) ) integral%char = char
                if( present(universe) ) integral%u = universe
                if( present(lps) ) integral%lps = lps
                if( present(halo) ) integral%halo = halo
                if( present(smith) ) integral%smith = smith
                if( present(galaxyBin) ) integral%galaxyBin = galaxyBin
                if( present(chi) ) integral%chi = chi

                integral%converged_2 = .false.

        end subroutine new_integral
        
        real(D) pure function getParam_integral(integral)
                type(tIntegral), intent(in) :: integral
                getParam_integral = integral%params(1)
        end function getParam_integral

        pure subroutine getParams_integral(integral, params)
                type(tIntegral), intent(in) :: integral
                real(D), intent(out) :: params(:)
                params = integral%params( 1 : size(params) )
        end subroutine getParams_integral

        pure subroutine getParams_2(integral, param1, param2)
                type(tIntegral), intent(in) :: integral
                real(D), intent(out) :: param1, param2
                param1 = integral%params(1)
                param2 = integral%params(2)
        end subroutine getParams_2

        pure subroutine getParams_3(integral, param1, param2, param3)
                type(tIntegral), intent(in) :: integral
                real(D), intent(out) :: param1, param2, param3
                param1 = integral%params(1)
                param2 = integral%params(2)
                param3 = integral%params(3)
        end subroutine getParams_3

        pure subroutine getParams_5(integral, param1, param2, param3, param4, param5)
                type(tIntegral), intent(in) :: integral
                real(D), intent(out) :: param1, param2, param3, param4, param5
                param1 = integral%params(1)
                param2 = integral%params(2)
                param3 = integral%params(3)
                param4 = integral%params(4)
                param5 = integral%params(5)
        end subroutine getParams_5

        pure subroutine getIntegerParams_2(integral, param1, param2)
                type(tIntegral), intent(in) :: integral
                integer, intent(out) :: param1, param2
                param1 = integral%integerParams(1)
                param2 = integral%integerParams(2)
        end subroutine getIntegerParams_2

        pure subroutine getIntegerParams_3(integral, param1, param2, param3)
                type(tIntegral), intent(in) :: integral
                integer, intent(out) :: param1, param2, param3
                param1 = integral%integerParams(1)
                param2 = integral%integerParams(2)
                param3 = integral%integerParams(3)
        end subroutine getIntegerParams_3

        character(len = INTEGRAL_CHAR_LEN) pure function getChar(integral)
                type(tIntegral), intent(in) :: integral
                getChar = integral%char
        end function getChar

        type(tHalo) pure function getHalo(integral)
                type(tIntegral), intent(in) :: integral
                getHalo = integral%halo
        end function getHalo

        type(tUniverse) pure function getUniverse_integral(integral)
                type(tIntegral), intent(in) :: integral
                getUniverse_integral = integral%u
        end function getUniverse_integral

        type(tLPS) pure function getLPS_integral(integral)
                type(tIntegral), intent(in) :: integral
                getLPS_integral = integral%lps
        end function getLPS_integral

        type(tGalaxyBin) pure function getGalaxyBin_integral(integral)
                type(tIntegral), intent(in) :: integral
                getGalaxyBin_integral = integral%galaxyBin
        end function getGalaxyBin_integral

        type(tSmith) pure function getSmith_integral(integral)
                type(tIntegral), intent(in) :: integral
                getSmith_integral = integral%smith
        end function getSmith_integral

        type(tChi) pure function getChi_integral(integral)
                type(tIntegral), intent(in) :: integral
                getChi_integral = integral%chi
        end function getChi_integral

        logical pure function converged_integral(integral)

                type(tIntegral), intent(in) :: integral

                integer :: i, j

                if(integral%dimension == 0) then
                        converged_integral = integral%converged
                else
                        converged_integral = .true.
                        do i = 1, integral%dimension
                                do j = 1, integral%dimension
                                        converged_integral = ( converged_integral .and. integral%converged_2(i, j) )
                                end do
                        end do
                end if

        end function converged_integral

        subroutine checkConvergence_integral(integral, errorMessage)
                type(tIntegral), intent(in) :: integral
                character(len = *), intent(in) :: errorMessage
                if( .not. converged(integral) ) then
                        write(0, *) errorMessage
                end if
        end subroutine checkConvergence_integral

        function converged_2(integral)
                type(tIntegral), intent(in) :: integral
                logical :: converged_2(integral%dimension, integral%dimension)
                converged_2 = integral%converged_2(1:integral%dimension, 1:integral%dimension)
        end function converged_2

        recursive subroutine integrate_0(integral, integrand, ans)

                use mPolynomialInterpolation

                type(tIntegral), intent(inout) :: integral
                interface
                        function integrand(integral, x)
                                use mPrecision
                                use mObject
                                real(D) :: integrand
                                type(tIntegral), intent(in) :: integral
                                real(D), intent(in) :: x
                        end function integrand
                end interface
                real(D), intent(out) :: ans
                real(D), parameter :: ONE_NINTH = 1.0_D / 9.0_D

                integer, parameter :: JMAXCLOSED = 20, & ! Max # of refinements for closed integration (20)
                                      JMAXOPEN = 14, &   ! Max # of refinements (14) :
                                      !JMIN = 5, &        ! Min # of refinements (5) :
                                      JMIN = 7, & ! MJ (5 was sometimes too small)
                                      ORDER = JMIN - 1 ! # of past refinements to include in interpolation (JMIN - 1):

                real(D), dimension( max(JMAXCLOSED, JMAXOPEN) + 1) :: h, s
                real(D) :: err, weight
                integer :: maxIters, j
!                real(D) :: ansArray(5), errArray(5)

                if(integral%open) then
                        maxIters = JMAXOPEN
                        weight = ONE_NINTH
                else
                        maxIters = JMAXCLOSED
                        weight = 0.25_D
                end if

                !print *,'integrate0: eps = ',integral%eps

                h(1) = 1.0_D
                iter_open: do j = 1, maxIters
                        if(integral%open) then
                                call midpoint( integral, integrand, j, s(j) )
                        else
                                call trapezoid( integral, integrand, j, s(j) )
                        end if
                        !print *,'j = ',j,'  s = ',s(j)
                        interp: if(j >= JMIN) then
                                call extrapolate( h(j - ORDER : j),  s(j - ORDER : j), ans, err)
                                !print *,'extrap: ans = ',ans,', err = ',err
!                                write(0, *) j, abs(err), '<=', integral%eps * abs(ans)
                                convergence: if( abs(err) <= integral%eps * abs(ans) ) then
                                        integral%converged= .true.
                                        return
                                end if convergence
                        end if interp
                        s(j + 1) = s(j)
                        h(j + 1) = h(j) * weight
                end do iter_open
                ! Some answer will be returned at this point, but it's hard to
                ! say how accurate it will be.
                integral%converged = .false.
        end subroutine integrate_0

        recursive subroutine integrate_2(integral, integrand, ans)

                use mPolynomialInterpolation
                use mMisc

                type(tIntegral), intent(inout) :: integral
                interface
                        subroutine integrand(integral, x, ans)
                                use mPrecision
                                use mObject
                                type(tIntegral), intent(in) :: integral
                                real(D), intent(in) :: x
                                real(D), intent(out) :: ans(:, :)
                        end subroutine integrand
                end interface
                real(D), intent(out) :: ans(:, :)

                real(D), parameter :: ONE_NINTH = 1.0_D / 9.0_D

                integer, parameter :: JMAXCLOSED = 20, & ! Max # of refinements for closed integration (20)
                                      JMAXOPEN = 14, &   ! Max # of refinements (14)
                                      JMIN = 5, &        ! Min # of refinements (5)
                                      ORDER = JMIN - 1   ! # of past refinements to include in interpolation (JMIN - 1)

                real(D), dimension( max(JMAXCLOSED, JMAXOPEN) + 1, size(ans, 1), size(ans, 2) ) :: h, s
                real(D) :: err( size(ans, 1), size(ans, 2) ), weight
                integer :: maxIters, j, m, n
!                real(D) :: ansArray(5), errArray(5)

                call assert( size(ans, 1) == size(ans, 2), &
                    'You have to generalize the integrate_2 subroutine (slightly) to allow for non-square integrands' )

                integral%dimension = size(ans, 1)

                if(integral%open) then
                        maxIters = JMAXOPEN
                        weight = ONE_NINTH
                else
                        maxIters = JMAXCLOSED
                        weight = 0.25_D
                end if

                h(1, :, :) = 1.0_D
                iter_open: do j = 1, maxIters
!                        write(0, *) j
                        if(integral%open) then
                                ! DEREK:
                                write(0, *) 'ERROR: Teach me to do open Romberg matrix integration.'
                                stop
!                                call midpoint( integral, integrand, j, s(j) )
                        else
                                call trapezoid_2( integral, integrand, j, s(j, :, :) )
                        end if
                        interp: if(j >= JMIN) then
                                do m = 1, integral%dimension
                                        inner: do n = 1, integral%dimension
                                                if( .not. integral%converged_2(m, n) ) then
                                                        call extrapolate( h(j - ORDER : j, m, n),  s(j - ORDER : j, m, n), ans(m, n), err(m, n) )
!            write(0, '(3i4, 2e11.3)') j, m, n, err(m, n), integral%eps * abs( ans(m, n) )
                                                        convergence: if(abs( err(m, n) ) <= integral%eps * abs( ans(m, n) )) then
                                                                integral%converged_2(m, n) =  .true.
!                                                                write(0, *) m, n, 'converged'
!                                                                write(0, '(11l3)') converged_2(integral)
!                                                                write(0, *)
                                                                if( converged(integral) ) return
!                                                                ans(m, n) = s(j, m, n)
                                                                cycle inner
                                                        end if convergence
!                                                else
!                                                        write(0, *) 'skipping', m, n
                                                end if
                                                s(j + 1, m, n) = s(j, m, n)
                                                h(j + 1, m, n) = h(j, m, n) * weight
                                        end do inner
                                end do
                        else interp
                                s(j + 1, :, :) = s(j, :, :)
                                h(j + 1, :, :) = h(j, :, :) * weight
                        end if interp
!                        s(j + 1, :, :) = s(j, :, :)
!                        h(j + 1, :, :) = h(j, :, :) * weight
                end do iter_open

                ! Some answer will be returned at this point, but it's hard to say how accurate it will be:
                ! DEREK
                integral%converged = .false. ! DEREK?
         end subroutine integrate_2

        subroutine quadrature(integrand, ll_ans, ul, x, w)

                interface
                        function integrand(x)
                                use mPrecision
                                real(D) :: integrand
                                real(D), intent(in) :: x
                        end function Integrand
                end interface
                real(D), intent(inout) :: ll_ans
                real(D), intent(in) :: ul
                real(D), intent(in) :: x(:), w(:)
                integer :: j
                real(d) :: dx, xm, xr
                xm = 0.5_D * (ul + ll_ans)
                xr = 0.5_D * (ul - ll_ans)
                ll_ans = 0.0_D
                do j = 1, size(w)
                        dx = xr * x(j)
                        ll_ans = ll_ans + w(j) * integrand(xm + dx)
                end do
                ll_ans = xr * ll_ans
        end subroutine quadrature
!C  (C) Copr. 1986-92 Numerical Recipes Software D041&0(9p#3.

        ! Upper and lower limits may be needed later,
        ! but now the subroutine Quadrature scales limits to -1,1
        pure subroutine gaussLaguerre(x,w)
                use mConstants
!LLUL            real(d), intent(in) :: ll,ul
                real(D), intent(out) :: x(:), w(:)
!PRECISION       real(D), parameter :: EPS = epsilon(1.0_D)
                real(D), parameter :: EPS = epsilon(1.0_R4)
                integer :: i, j, m, n
                real(D) :: p1, p2, p3, pp, xl, xm, z, z1
                n = size(x)
                m = (n + 1) / 2
!LLUL            xm=0.5_d*(ul+ll)
!LLUL            xl=0.5_d*(ul-ll)
                do i = 1, m
                        z = cos(PI * ( real(i, D) - 0.25_D ) &
                                 / ( real(n, D) + 0.5_D ))
                        do
                                p1 = 1.0_D
                                p2 = 0.0_D
                                do j = 1, n
                                        p3 = p2
                                        p2 = p1
                                        p1 = (( 2.0_D * real(j, D) - 1.0_D ) &
                                               * z * p2 &
                                                - ( real(j, D) - 1.0_D ) * p3) &
                                             / real(j, D)
                                end do
                                pp = n * (z * p1 - p2) / (z * z - 1.0_D)
                                z1 = z
                                z = z1 - p1 / pp
                                if( abs(z - z1) <= EPS ) exit
                        end do
!LLUL                    x(i) = xm - xl * z
                        x(i) = -z
!LLUL                    x(n + 1 - i) = xm + xl * z
                        x(n + 1 - i) = z
!LLUL                    w(i) = 2.0_D * xl / ( (1.0_D - z * z) * pp * pp )
                        w(i) = 2.0_D / ( (1.0_D - z * z) * pp * pp )
                        w(n + 1 - i) = w(i)
                end do
        end subroutine gaussLaguerre
!C  (C) Copr. 1986-92 Numerical Recipes Software D041&0(9p#3.

        !**********************************************************************
        ! This calculates the order-th stage of refinement using
        ! the extendend trapezoidal rule. It requires the integrand to be
        ! finite at the limits ll,ul. When called with order = 1, the crudest
        ! estimate of the integral is returned. Subsequent calls with
        ! order = 2, 3, ... (should be called in sequence) will improve the accuracy
        ! of the answer ("ans") by adding 2 ** (n - 2) additional interior points.
        !**********************************************************************
        recursive subroutine trapezoid(integral, integrand, order, ans)
                type(tIntegral), intent(in) :: integral
                interface
                        function integrand(integral, x)
                                use mPrecision
                                use mObject
                                real(D) :: integrand
                                type(tIntegral), intent(in) :: integral
                                real(D), intent(in) :: x
                        end function integrand
                end interface
                integer, intent(in) :: order
                real(D), intent(inout) :: ans
                integer :: it, j
                real(D) :: del, sum, oneOver_tnm, x, newLL, newUL

                call changeVariables(integral, newLL, newUL)

                crude: if(order == 1) then
                        crude_change: select case(integral%changeVariables)
                        case default crude_change ! no change of variables
                                ans = 0.5_D * (newUL - newLL) * ( integrand(integral, newLL) + integrand(integral, newUL) )
                        case('EXP') crude_change ! x -> e^x
                                ans = 0.5_D * (newUL - newLL) &
                                    * (exp(newLL) * integrand(integral, exp(newLL) ) + exp(newUL) * integrand( integral, exp(newUL) ))
                        end select crude_change
                else crude
                        it = 2 ** (order - 2)
                        oneOver_tnm = 1.0_D / it
                        del = (newUL - newLL) * oneOver_tnm
                        x = newLL + 0.5_D * del
                        sum = 0.0_D
                        change: select case(integral%changeVariables)
                        case default change
                                sum_none: do j = 1, it
                                        sum = sum + integrand(integral, x)
                                        x = x + del
                                end do sum_none
                        case('EXP') change ! x->e^x
                                sum_exp: do j = 1, it
                                        sum = sum + exp(x) * integrand( integral, exp(x) )
                                        x = x + del
                                end do sum_exp
                        end select change
                        ans = 0.5_D * ( ans + (newUL - newLL) * sum * oneOver_tnm )
                end if crude
        end subroutine trapezoid

        recursive subroutine trapezoid_2(integral, integrand, order, ans)
                type(tIntegral), intent(in) :: integral
                interface
                        subroutine integrand(integral, x, ans)
                                use mPrecision
                                use mObject
                                type(tIntegral), intent(in) :: integral
                                real(D), intent(in) :: x
                                real(D), intent(out) :: ans(:, :)
                        end subroutine integrand
                end interface
                integer, intent(in) :: order
                real(D), intent(inout) :: ans(:, :)
                integer :: it, j
                real(D) :: del, oneOver_tnm, x, newLL, newUL
                real(D), dimension( size(ans, 1), size(ans, 2) ) :: ansLL, ansUL, sum
!                real(D), dimension( 0:size(ans, 1) - 1, 0:size(ans, 1) - 1 ) :: ansLL, ansUL, sum

                call changeVariables(integral, newLL, newUL)

                crude: if(order == 1) then
                        crude_change: select case(integral%changeVariables)
                        case default crude_change ! no change of variables
!                                ans = 0.5_D * (newUL - newLL) * ( integrand(integral, newLL) + integrand(integral, newUL) )
                                call integrand(integral, newLL, ansLL)
                                call integrand(integral, newUL, ansUL)
                                ans = 0.5_D * (newUL - newLL) * ( ansLL + ansUL )
                        case('EXP') crude_change ! x -> e^x
!                                ans = 0.5_D * (newUL - newLL) &
!                                    * (exp(newLL) * integrand(integral, exp(newLL) ) + exp(newUL) * integrand( integral, exp(newUL) ))
                                call integrand( integral, exp(newLL), ansLL )
                                call integrand( integral, exp(newUL), ansUL )
                                ans = 0.5_D * (newUL - newLL) * ( exp(newLL) * ansLL + exp(newUL) * ansUL )
                        end select crude_change
                else crude
                        it = 2 ** (order - 2)
                        oneOver_tnm = 1.0_D / it
                        del = (newUL - newLL) * oneOver_tnm
                        x = newLL + 0.5_D * del
                        sum = 0.0_D
                        change: select case(integral%changeVariables)
                        case default change
                                sum_none: do j = 1, it
                                        call integrand(integral, x, ansLL)
                                        where( .not. integral%converged_2(1:size(ans, 1), 1:size(ans, 2)) )
                                                sum = sum + ansLL
                                        end where
                                        x = x + del
                                end do sum_none
                        case('EXP') change ! x->e^x
                                sum_exp: do j = 1, it
                                        call integrand( integral, exp(x), ansLL )
                                        where( .not. integral%converged_2(1:size(ans, 1), 1:size(ans, 2)) )
                                                sum = sum + exp(x) * ansLL
                                        end where
                                        x = x + del
                                end do sum_exp
                        end select change
                        ans = 0.5_D * ( ans + (newUL - newLL) * sum * oneOver_tnm )
                end if crude
        end subroutine trapezoid_2

        !***************************************************************************
        ! This calculates the order-th stage of refinement using
        ! the extendend midpoint rule (does not require integrand to be finite
        ! at the limits ll, ul. When called with order = 1, the crudest
        ! estimate of the integral is returned. Subsequent calls with
        ! order = 2, 3... (should be called in sequence) will improve the accuracy
        ! of the answer ("ans") by adding 2 * 3 ** n - 2 additional interior points.
        !***************************************************************************
        recursive subroutine midpoint(integral, integrand, order, ans)
                use mConstants
                type(tIntegral), intent(in) :: integral
                interface
                        function integrand(integral, x)
                                use mPrecision
                                use mObject
                                real(D) :: integrand
                                type(tIntegral), intent(in) :: integral
                                real(D), intent(in) :: x
                        end function integrand
                end interface
                integer, intent(in) :: order
                real(D), intent(inout) :: ans
                integer :: it, j
                real(D) :: ddel, del, sum, oneOver_tnm, x, newLL, newUL

                call changeVariables(integral, newLL, newUL)

                crude: if(order == 1) then
                        x = 0.5_D * (newUL + newLL)
                        crude_change: select case(integral%changeVariables)
                        case default crude_change ! no change
                                ans = (newUL - newLL) * integrand(integral, x)
                        case('INV') crude_change ! x->1/x
                                ans = (newUL - newLL) / (x * x) &
                                      * integrand(integral, 1.0_D / x)
                        case('EXP') crude_change ! x->e^x
                                ans = (newUL - newLL) * exp(x) &
                                      * integrand(integral, exp(x))
                        end select crude_change
                else crude
                        it = 3 ** (order - 2)
                        oneOver_tnm = 1.0_D / it
                        del = (newUL - newLL) * oneOver_tnm * ONE_THIRD
                        ddel = del + del
                        x = newLL + 0.5_D * del
                        sum = 0.0_D
                        change: select case(integral%changeVariables)
                        case default change ! no change
                                sum_none: do j = 1, it
                                        sum = sum + integrand(integral, x)
                                        x = x + ddel
                                        sum = sum + integrand(integral, x)
                                        x = x + del
                                end do sum_none
                        case('INV') change ! x->1/x
                                sum_inf: do j = 1, it
                                        sum = sum &
                                            + integrand(integral, 1.0_D / x)&
                                              / (x * x)
                                        x = x + ddel
                                        sum = sum &
                                            + integrand(integral, 1.0_D / x) &
                                              / (x * x)
                                        x = x + del
                                end do sum_inf
                        case('EXP') change ! x->e^x
                                sum_exp: do j = 1, it
                                        sum = sum + exp(x) &
                                            * integrand( integral, exp(x) )
                                        x = x + ddel
                                        sum = sum + exp(x) &
                                            * integrand( integral, exp(x) )
                                        x = x + del
                                end do sum_exp
                        end select change
                        ans = ( ans + (newUL - newLL) * sum * oneOver_tnm ) &
                              * ONE_THIRD
                end if crude
        end subroutine midpoint

        subroutine changeVariables(integral, changedLL, changedUL)
                type(tIntegral), intent(in) :: integral
                real(D), intent(out) :: changedLL, changedUL
                
                select case(integral%changeVariables)
                case default ! no change---do nothing
                        changedLL = integral%ll
                        changedUL = integral%ul
                case('INV') ! x->1/x

                        if(.not. integral%open) then
                                write(0, *) 'ERROR: Currently I can only do'
                                write(0, *) '       x->1/x change of variables'
                                write(0, *) '       with open integration'
                                write(0, *) '       method.'
                                stop
                        end if

                        changedLL = 1.0_D / integral%ul
                        changedUL = 1.0_D / integral%ll
                case('EXP') ! x->e^x, OK

                        if(integral%ll == 0.0_D) then
                                write(0, *) 'ERROR: Lower limit should not be'
                                write(0, *) '       zero when making change'
                                write(0, *) '       of variables x->e^x.'
                                stop
                        end if
                        if(integral%ul == 0.0_D) then
                                write(0, *) 'ERROR: Upper limit of integral'
                                write(0, *) '       cannot be zero when making'
                                write(0, *) '       making change of variables'
                                write(0, *) '       x->e^x'
                                stop
                        end if

                        changedLL = log(integral%ll)
                        changedUL = log(integral%ul)
                end select
        end subroutine changeVariables

end module mIntegration
